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ABSTRACT 

We develop a novel method of measuring the lensing distortion profiles of clusters with stacking the “scaled” 
amplitudes of background galaxy ellipticities as a function of the “scaled” centric radius according to the 
Navarro-Frenk-White (NFW) prediction of each cluster, based on the assumption that the different clusters in 
a sample follow the universal NFW profile. First we demonstrate the feasibility of this method using both the 
analytical NFW model and simulated halos in a suite of high-resolution /V-body simulations. We then apply, 
as a proof of concept, this method to the Subaru weak lensing data and the XMM/Chandra X-ray observables 
for a sample of 50 massive clusters in the redshift range 0.15 < z < 0.3, where their halo masses differ from 
each other by up to a factor of 10. To estimate the NFW parameters of each cluster, we use the halo mass 
proxy relation of X-ray observables, based on either the hydrostatic equilibrium or the gas mass, and then infer 
the halo concentration from the model scaling relation of halo concentration with halo mass. We evaluate a 
performance of the NFW scaling analysis by measuring the scatters of 50 cluster lensing profiles relative to 
the NFW predictions over a range of radii, 0.14 < R/[h 'Mpc] < 2.8. We found a 4 - 6<x level evidence of 
the universal NFW profile in 50 clusters, for both the X-ray halo mass proxy relations, although the gas mass 
appears to be a better proxy of the underlying true mass. By comparing the measurements with the simulations 
of cluster lensing profiles taking into account the statistical errors of intrinsic galaxy shapes in the Subaru data, 
we argue that additional halo mass errors or intrinsic scatters of cr( M ^qo c ) / M ^oo c ~ 0.2 - 0.3 could reconcile a 
difference between the measurements and the simulations. This method allows us to some extent to preserve 
characteristics of individual clusters in the statistical weak lensing analysis, thereby yielding a new means of 
exploiting the underlying genuine form of the halo mass profile and the halo mass proxy relations via weak 
lensing information, under the assumption of the existence of the universal profile. 

Subject headings: cosmology: observations - dark matter - galaxies: clusters: general - gravitational lensing: 
weak 


1. INTRODUCTION 

Clusters of galaxies are the largest, gravitationally bound 
objects in the Universe, and the formation and evolution pro¬ 
cesses are dominated by gravitational effects mainly due to 
dark matter. Hence clusters provide us with a useful labora¬ 
tory for studying the nature of dark matter (Clowe et al.|2006 1 
as well as constraining cosmology, e.g. from the abundance 
of clusters found from a survey volume (Vikhlinin et al. 2009 


|Rozo et al.|2010| |Oguri & Takada|2011) . However, to attain 

the full potential of cluster-based cosmology from upcoming 
wide-area surveys such as the Subaru Hyper Suprime-Canrl 
and the Dark Energy Surve\0 requires a further understand^ 
ing of the physical processes in clusters. 

One of the most important predictions in /V-body simula¬ 
tions of the A-dominated, cold dark matter structure forma¬ 
tion model (ACDM) is the emergence of universal mass den¬ 
sity profile - that is, the mass density profile of dark matter 
halos can be well fitted by a “universal” two-parameter fam¬ 
ily of the model profile over a wide range of halo masses, 
first proposed in Navarro et aTT| ( 1996| 1997 hereafter NFW). 
The NFW profile predicts a monotonically steepened profile 
with increasing radius, with logarithmic slopes shallower than 


1 http://www.naoj.org/Projects/HSC/j_index.html 

2 http://www.darkenergysurvey.org 


an isothermal sphere interior to the characteristic “scale” ra¬ 
dius r < r s , but steeper at larger radius, approaching r 3 at 
the virial radius, r —» r VH (see also |Dalal et al.|2010| for dis¬ 
cussion on the physical origin within the framework of the 
hierarchical ACDM model). Further, the ratio of the charac¬ 
teristic scale radius to the virial radius, which characterizes 
the degrees of central concentration of the mass distribution 
- the so-called halo concentration c - tends to be lower for 
more massive halos. In addition the halo concentration of a 
fixed halo mass displays intrinsic scatters typically given by 
CT] nc ~ 0.2, originating from details of the mass accretion or 


mation (Bullock et al. 20011 Wechsler et al.|2002; 

Zhao et al. 

2003 

Duffy et al. 2008 Zhao et al. 2009; Bhattac 

tarya et al. 


dark matter halos are important predictions of ACDM model, 
and need to be carefully tested by measurements. 

Gravitational lensing is a unique, powerful method enabling 
one to probe the matter distribution in galaxy clusters irre¬ 
spective of their physical and dynamical states (Schneider] 


2006f. Several works have investigated the mass density pro 


file over a wide range of radii by combining the strong and 
weak lensing at the small and large radii, observationally the 
exquisite high-resolution images of Hubble Space Telescope 
and the wide-held ground based telescopes such as the Subaru 
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Telescope (Tyson & Fischer 

1995 

Cneib et al. 

2003 Broad- 

hurst et al. 2005 1 Smith et"al 

. 2005 

Oguri et"al 

. 2005, 2012, 


Merten et al. 2014 1 . In addition, the stacked weak lensing 

analysis combining a sample of clusters has been proven to 
be a robust, powerful method of probing the average mass 
distribution of the sampled cluster s (Johnston et al. 2007 Ok¬ 
ie et al.|2010a[ Oguri et al.|20121|O 
et al.|2014b. These works have s 


kabe et al.|2013||tJmetsu 

own that the average mass 

profile measured from the stacked lensing is in remarkably 
nice agreement with the NFW prediction. Another advantage 
of stacked lensing is it allows one to probe the mass distribu¬ 
tion even for less massive halos, such as galaxy-scale halos, 
as long as a sufficient number of sam pled halos (e.g. galaxies) 
are used in the analy sis (Mandelbaum et al.|2005| Leauthaud 


|et al.|2010l|Miya take et a l.|2013[ ). However, a downside of the 

stacked lensing method is the loss of lensing information of 
individual clusters. Hence a knowledge of the distribution of 
the underlying halo parameters in the sampled clusters such as 
their halo masses is of critical importance in order not to have 
any bias in the NFW parameters inferred from the stacked 
lensing signals ( |Oguri & Takada|2011[> . This is equivalent to 
the importance of exploring a well-calibrated proxy relation 
of cluster observables with halo mass or more generally halo 


parameters (Rozo et all2009 

Zhang et al.|2010 ; Okabe et al. 

2010 

Zhang et al. 20 IT 

; Ma 

hdavi et al. 2013 

von der Linden 

et al. 

2014, Martino et a 

. 2014; Okabe et al. 

2014; Hoekstra 

et al. 

2015). 


The purpose of this paper is to develop a novel method 
of measuring the lensing distortion profiles of clusters, mo¬ 
tivated by the NFW prediction. We propose the “NFW scal¬ 
ing” analysis for the lensing measurements, which is done by 
averaging the “scaled” amplitudes of background galaxy el- 
lipticities in each bin of the “scaled” radii according to the 
NFW prediction of individual cluster. With this NFW scaling 
method, we can address whether clusters display the univer¬ 
sality of their lensing profiles as seen in simulations. First, 
to demonstrate the feasibility of the NFW scaling analysis, 
we will use simulations of cluster lensing observables based 
on a suite of high-resolution /V-body simulations. Then, as 
a proof of concept, we will apply this method to a sample 
of 50 massive clusters in the redshift range 0.15 < z < 0.3 
that have been observed with the Subaru telescope by the Lo- 


cal Cluster Substructure Survey (LoCu 

SS, and also see Okabe 
for details!]] Note that 

et al.]2010a 

2013 

Martino et al.12014 

this study is basec 

on the published results of LoCuSS, and is 


not performed within the collaboration. To estimate the NFW 
scaling of each cluster, we will use the halo mass estimate in 
|Martino et al.| ( ]2014} > based on the XMM and/or Chandra X- 
ray observables, and use the halo concentration inferred from 
the model scaling relation between halo mass and concentra¬ 
tion in Diemer & Kravtsov] ( )2015) >. Then by comparing the 
scatters of 50 cluster lensing profiles relative to the NFW pre¬ 
dictions for two cases with and without the NFW scaling, we 
test the performance of this method as well as the universality 
of the cluster mass distribution. 

The structure of this paper is as follows. In §|2] after briefly 
reviewing the lensing observables of NFW halo] we will de¬ 
rive an estimator of the lensing distortion profile measure¬ 
ment with NFW scaling. Then we study the feasibility of this 

’ Based in part on data collected at Subaru Telescope and obtained from 
the SMOKA, which is operated by the Astronomy Data Center, National As¬ 
tronomical Observatory of Japan. 


method using analytical NFW models and /V-body simula¬ 
tions. In §[3] we first describe the Subaru weak lensing catalog 
and the X-ray observables for the sample of massive clusters 
we use in this paper, and show the results of the application 
of this method to the Subaru data. § |4] is devoted to discus¬ 
sion and conclusion. Unless stated otherwise, we will adopt a 
flat ACDM cosmology with O,,, = 0.27, Oa = 0.73, and the 
Hubble parameter h = Hq/( 100 km s l Mpc ') = 0.70. 

2. METHODOLOGY: STACKED WEAK LENSING WITH 
NFW SCALING 

2.1. Lensing ofNavarro-Wliite-Frenkhalo 

The Navarro-Frenk-White (1997; hereafter NFW) mass 
density profile for a halo is parametrized by two parameters 
as 

Pc 


Pnfw(t) 


(r/r. s )(l + r/r s ) 2 ' 


(1) 


where r s is the scale radius and p c is the central density param¬ 
eter. The parameter p c is specified by imposing that the mass 
enclosed within a sphere of a given overdensity A is equal to 
the halo mass M A , 


Ap cr (z)c 




3ot N fw(c a ) 4 ^TjW NFW (ca) 


( 2 ) 


where ot nfw (ca) = Jf A dx x/( 1 +x ) 2 = ln(l +c A )-c A /(l + c A ), 
c A = ta/Cs, a concentration parameter, and A(z) is a nonlin¬ 
ear overdensity introduced to define the interior mass for each 
halo. Note that throughout this paper we employ halo mass 
definition with respect to the critical density, not the mean 
mass density: = (47r/3)r^p cr (z)A. Although we focus on 

the NFW profile in this paper, the method developed in this 
paper can be applied to any other univeral profile such as gen¬ 
eralized NFW or Einasto profile (Einasto 1965l|Merritt et al. 

rot . 

Several works have shown a scaling relation of the halo 
concentration with halo mass, using numerical simulations or 
based on analytical arguments (|Bullock et al.|2001| [WechsTer 

38 ) i 


et al.l2002tlZhao et al.12003][Puff y et al. 2008) (most recently 

Diemer & Kravtsov 2015 hereafter DK15 and see references 


therein). As for our fiducial model, we adopt the publicly- 
available code provided by B. Diemer to compute the halo 
mass and concentration relation in DK15. Note that we used 
the “median” relation, rather than the mean, for our default 
choice as recommended in DK15. The mass estimates from 
the X-ray observables are not M 2000 and rather the interior 
mass of a greater overdensity such as M 5 () o c . Assuming that 
a halo exactly follows the NFW profile, we can convert the 
scaling relation calibrated for M 200 C to the C 5 qo c -M 5 oq c rela¬ 
tion, based on the method in |Hu & Kravtsov| (|2003). The 
public code of DK15 allows us to compute the halo concen¬ 
tration for an input overdensity based on this method. We will 
also study how possible variations in the c-M relation affect 
the results of this paper. 

For an NFW profile, we can derive an analytical exp ressi on 
for the lensing convergence and shear profiles (Bartelmann 
1996| [Golse"& Kneib|2002T >: 


k nfw (/?) = 




(R) 


Ecr(Zi, Zs) 

A Z NFW (R) 
£ cr (Zl,Z s ) 


2 p c r. 


= 2p c r s 


g NFW (R/r s ) 

^-crit(Zl, Zs) 

/ NFW (R/r s ) 


^-crit(Zl, Zs) 


(3) 
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where R is the projected comoving radius from halo center, 
and the functions / NFW (x) and g NFW (x) are given by 


1- cosh (x < 1) 

■2 X 


g NFW (x) = 


x 2 -i\ VT^x- 


* - 1 \ Vx 2 - 1 


1 -tl 

iCOS - 


and 


/ NFW (x) 


1 


2 x 

x 2 n 2 + 1 — r 2 
- -21n2, 


2 x 

— In - 


1 


x 2 - 1 


(x = 1) (4) 

(x> 1) 


, 2 - 3x 2 , 1 \ „ 

1 +- cosh - , (x < 1) 

X 2 V 1 - X 2 x ) 

(x = 1)(5) 

, 2 - 3x 2 _! 1 \ 

1 +- cos - , (x > 1). 

x 2 Vx 2 — 1 X ) 


The critical surface mass density X cr ; t for a given system of 
lens cluster and source at redshifts z\ and z s , respectively, is 
given as 


2 iciit(Zl, Zs) — 


Da(Zs) 


4 nG D a (z\)D a (zi, Zs)( 1 + zi) 2 ’ 


( 6 ) 


where D A (z) is the angular diameter distance and the factor 
(1 + si) 2 is from our use of the comoving scale. From Eqs. 0 
and Q, we can find that the lensing amplitudes of an NFW 
halo scale with the NFW parameters (M A , c A ) as 


^fw^fw cc 2p cTs oc M a /(r 2 m NFW (c A )) 


ocM/ c 2 A /m NFW (c A ). 


(7) 


If we employ the c A -M A scaling relation given as c A (M A ) oc 
M “, the lensing amplitudes roughly scale with halo mass as 
y+ FW p,- t/3-2 a as f unc tion m(c A ) has a weak dependence 
on halo mass. Note that, since the cluster sample is among the 
most massive clusters, we have checked that the 2-halo term 
is much smaller than the above 1-halo term, by a factor of 
100, over a range of the radii we consider (e.g., see Oguri & 
|Takada|201 1| Takada & Spergel|2014| i. Therefore we ignore 
the 2-halo term for the following analysis. 

An actual lensing observable estimated from ellipticities of 
background galaxies for an NFW lens is the lensing “distor¬ 
tion” profile or reduced shear profile: 


<«+>(*)■ 


1 -k nfw <7?)’ 


( 8 ) 


where e + is the tangential component of the ellipticities with 
respect to cluster center. The reduced shear correction is not 
negligible at the inner radii, and we need to take into account 
the correction. 


2.2. Stacked lensing without NFW scaling 

For the standard method to estimate the stacked lensing pro- 
filefor N c clusters, we follow the m ethod in Johnsto n et al.| 
( |2007| ) and |Mandelbaum et al.| ( |2013l >: 

_ j Nc 

<A2 )(R) = - z z 

t s a .IU',71 V. kt? 


where e (Vi)+ is the tangential ellipticity of the s fl -th back¬ 
ground galaxy in the o-th cluster region, and N is the nor¬ 
malization factor defined as 


N c 

N = TjYj W ^> (10) 

< 2=1 S a 

The summation £„ runs over the sampled clusters, from a - 1 
to N c , and the summation Zs„:|r, |es mns over all the back¬ 
ground galaxies that reside in the annulus of radius R from 
the a-th cluster center to within the bin width. We employ the 
weight given as 


W(a,* 0 ) 


1 

Zcr(Za,zJ 2 (el a) + cr 2 (Sa)e + a 2 )' 


( 11 ) 


where z Sa is the redshift of the ,v n -th background galaxy, e (Vi) is 
the ellipticity amplitude, o~( Sa)e is the measurement error and 
a is the constant factor to regularize the weight for which we 
adopt a = 0.4 (|Okabe et al.||2010a| . Note that we employ 
the average redshift for all the source galaxies in each cluster 
region, as described below in detail. 

Since we need to employ a finite number of the radial bins 
to study the “shape” of lensing distortion profile, which bin¬ 
ning scheme to use is not so clear. As for the representative 
value of a given radial bin, we use the average of radii of back¬ 
ground galaxies that reside in the annulus taking into account 
their weights: 


2 a =l 2s 0 ;|R (a)j j6« W(a,s a )R(a)s a 

K ^-• (TA) 

In the literature the area-weighted value of each radial bin is 
often used. We have checked that, using an analytical NFW 
profile and taking the actual distribution of background galax¬ 
ies in the Subaru data, the above radial binning is more accu¬ 
rate in the sense that the distortion profile is in better agree¬ 
ment with the model NFW profile amplitude inferred by the 
representative value of the radial bin, less than 1% in the frac¬ 
tional difference for most cases. 

The statistical uncertainty of the stacked lensing at each ra¬ 
dial bin can be estimated as 

j N c 

0-(A #) 2 ^X 2 W i,s a )^r( a )(Za,Zsfel a y d3) 

a=l \cR 


In this paper we consider the intrinsic ellipticities as a source 
of the statistical errors in the lensing measurement, and ignore 
the cosmic shear contribution that arises from different mass 
distribution along the same line of sight to the cluster. For 
the application of this method to the Subaru data that we will 
show below, this is a good approximation in practice, because 
the mean number density of background galaxies is small, 
about 5 arcmin 2 , after a secure s electio n of background “red” 
galax ies as we will discuss in § |3.1.1 (also see Okabe et al. 
|2013| for the details). 

When comparing the measured lensing profile to an NFW 
model, we need to account for the contribution of reduce 
shear. In this paper, assuming that all the clusters follow a sin¬ 
gle NFW profile in average sense, we model the the stacked 












































4 


Niikura et al. 


lensing profile, according to Eqs. ([3]) and ([8]l as 


(AS) (R) 


AS nfw ( R ) 

1 - k nfw (R) 


AS nfw ( R ) 


1 + ( — ) S NFW «fl» 


, (14) 


where the notation “ <=> ” is meant to denote the compari¬ 
son between the measurement (left-hand side) and the model 
profile (right-hand side). The notation ( ) on the right-hand 

side denotes the average taking into account the weights of 
background galaxies in each cluster region as in Eq. ( [12| ). We 
will use the above equation to estimate the halo mass ana con¬ 
centration parameter, M A and c A , from the measured lensing 
profile. 


2.3. Stacked lensing with NFW scaling 

Now we consider the stacked lensing analysis with “NFW 
scaling”. To implement this method we combine the weak 
lensing measurement and X-ray observables, where the X-ray 
observables are needed to estimate halo mass of each clus¬ 
ter independently of the lensing observables. Assuming that 
each of the sampled clusters follows an NFW profile specified 
by their respective parameters, M( fl > and C( a) , we can define 
an estimator of the normalized NFW lensing profile from the 
measured ellipticities of background galaxies, as motivated by 

Eq® 


</»»>(*) = ^ Z Z 

a=l s a ;|x (a)ifl |e* 


W(a, ) 


(15) 

Here M, x , and c:, x , are the halo mass and concentration for 
the a-th cluster, estimated from the X-ray observables (see 
below for details). The scaled radius in the above equation, 
x, is defined for the a-cluster as = R( a)So / r s (M * a) , c x () ), 
where r s is the scale radius of NFW profile, r s = r A /c A . We 
use the representative value of each radial bin, x, defined in 
a similar manner to Eq. ( fT2] >. The central density parameter 
of NFW profile, p c , can be estimated from ;V/ X . and c* a) for 

the a-th cluster, from Eq. ({2J1. Note that the profile (/ NFW ) 
and the radius x are dimension-less. With the above NFW 
scaling, weak lensing signals due to less massive halos than 
the mean mass in the sampled clusters are up-weighted, while 
the signals of more massive halos are down-weighted. 

Similarly, the measurement errors of the stacked profile at 
each radial bin are estimated as 


CT ^yNFW^(X') 


1 


N c 

2 N 2 ^f 


Z 


W M T. a (a)(Za,Zj e (sJ 


a=l 4 P c{Mf ay Cfj cf fl) )' 


( 16 ) 

To test an improvement in the stacked lensing analysis of 
NFW scaling compared to the standard stacked lensing, we 
compare the scatters of lensing distortion profiles of the sam¬ 
pled clusters relative to the NFW prediction. To be more pre¬ 
cise, we quantify the scatters by 


" ! = ZZ 


a= 1 i 


\NZ (a) (R {a)l ) - AS bf - NFW (fl (fl)i ; M bf , c bf )] 2 


(17) 


or 


N c 

^w-scaling — ^ j ^ ^ 
ci— 1 i 




0-fNFW (a) (X( a)i ) 2 


(18) 


Here AS (aJ and / XFW are the measured distortion profile with¬ 
out and with NFW scaling for the a-th cluster, which are es¬ 
timated in the similar manner to Eqs. d9| and dl5l, and (r A v (fl) 
and (r^M-w (H) are the errors at each radial bin, estimated simi¬ 
larly to Eqs. (13 1 and ( |l6| , respectively. AS hf ~ NFW (f?) is the 
best-fit NFW profile of the stacked lensing profile (Eq. p~4] ). 
For the NFW scaling case, we similarly include the reduced 
shear correction: we multiply the function f NFW (x) (Eq.pji by 


the function, 1 + (1/E cr ) S NFW (x) as in Eq. (14), where 


we used the best-fit NFW model of the stacked lensing pro¬ 
file without NFW scaling in order to compute E NFW (x). The 
above d 2 and d 2 ; sca | in „ are equivalent to the log-likelihood 
functions of lensing distortion profiles of N c clusters assuming 
that the statistical errors are given by the intrinsic ellipticities. 
The radial bin R {ll ]l or x (a)! for the a-th cluster is similarly com¬ 
puted by Eq. ( p~2] > from the background galaxies that reside in 
the annulus ofthe cluster. However, comparing the scatters 
of lensing distortion profiles with and without NFW scaling 
requires a careful treatment of the radial binning. As we will 
later describe in more detail, we will employ the radial bin¬ 
ning scheme so as to preserve the same background galaxies 
in the ;-th radial bins with or without the NFW scaling. With 
this binning scheme, the relation d 2 = ^_ sca | in „ holds if set¬ 
ting the model profiles to A£ bf ~ NFW = / NFW = 0 Q If the 
lensing distortion profiles of the sampled clusters are similar 
in their shapes and amplitudes, following the NFW profile. 


the value of d~ should become smaller: d 


w-scaling 


<d 2 . 


2.4. Testing the method with N-body simulations 

In this subsection, before going to the Subam data, we test 
our method using analytical NFW model and high-resolution 
A'-body simulations. For the sake of convenience to compare 
with the following sections, we consider 50 clusters in this 
section as the 50 Subaru clusters. 

First let’s consider an ideal case, albeit unrealistic, that 
each of 50 clusters exactly follows an NFW profile. Fig¬ 
ure [I] shows the lensing profiles with or without the NFW 
scaling for 50 halos. To take into account variations in halo 
masses that resemble the 50 clusters, we assign one-by-one 
the X-ray inferred masses of 50 cluster to NFW halo^] Note 
that we use the c-M scaling relation in DK15 to compute the 
halo concentration for each NFW halo. The different blue 
curves show each NFW distortion profile relative to the best- 
fit NFW model of the stacked distortion profile, as a func¬ 
tion of the radius relative to the scale radius of the best-fit 
NFW model. Here we consider the same range of radii, 
0.14 < f?/[/i _1 Mpc] < 2.8 for all the halos as we will do 
for actual analysis of Subaru data. For the range of clus¬ 
ter masses, the lensing distortion amplitudes differ from each 
other by up to a factor of 5. 

4 If we set A£ bf - NFW = y NFW = 0 in Eqs. j i l\ and j 1 8} , the values of d 2 
give the cumulative signal-to-noise ratio of the lensing distortion measure¬ 
ments for the N c clusters. 

5 W e here employed the hydrostatic equilibrium mass in |Martino et al.| 
|2Q14| , which was estimated from the X-ray observables of each cluster. 
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Figure 1 . The distribution of NFW lensing profiles for 50 halos for each of 
which we took the X-ray inferred mass of 50 Subaru clusters (here the hy¬ 
drostatic equilibrium mass in Table QJ and assumed the halo concentration 
based on the halo m ass and concentration relation, c = c(Msoo c ), infDiemer 
& Kravtsov (2015 hereafter DK15). The blue curves are the lensing profiles 
without “NFW scaling”, i.e. the standard method, but each curve is normal¬ 
ized by the best-fit NFW profile to the stacked profile of 50 profiles, and is 
plotted as a function of the radius relative to the scale radius of the best-fit 
NFW model. The red curves are the lensing profiles with “NFW scaling”, 
computed assuming that the halo mass and concentration of each halo are a 
priori known — an ideal case. Note that we fixed the same range of radii, 
0.14 < /?/[/* _1 Mpc] < 2.8, in the comoving length units for both the results. 


On the other hand, the red curves in Figure [I] show the pro¬ 
files after the NFW scaling implementation, assuming that 
the true mass and concentration of each cluster are a priori 
known, i.e. an ideal case. Each curve is the fractional pro¬ 
file relative to the NFW distortion profile including the re¬ 
duced shear correction, / NFW (x) (Eq.]5]and see below Eq. 15 i. 
The deviation from unity is due to an imperfect correction 
of the reduced shear: the nonlinear correction becomes non- 
negligible at small radii, and breaks the universality of the 
NFW lensing profile. The horizontal axis is in the units of 
the “scaled” radius, R/r s , where r s is the NFW scale radius of 
each halo. Due to the radial transformation from the original 
fixed range of R, the range of the scaled radius x, covered by 
each halo, differ from each other. The figure shows that the 
NFW scaling significantly reduces the scatters of lensing pro¬ 
files, making the differences within 20% over a range of radii 
we consider. 

Obviously actual clusters have much more complicated 
mass distribution than an analytical NFW model: intrinsic 
scatters of halo concentration, aspherical mass distribution, 
substructures and so on. To study these effects we use simu¬ 
lated halos of cluster scales, generated from a high-resolution 
A'-body simulation in Takahashi et al.|(|2012|). In brief the N- 
body simulation was ran with the publicly-available Gadget- 


2 code (|Springel et a1. 11 2001 |Springel 2005) assuming the 
WMAP cosmology. The simulation employed 1024 3 particles 
in a box of 320 h 1 Mpc on a side. The mass resolution (the 
particle mass) is 2.3 x 10 9 lr x M & , so is sufficient to resolve 
cluster-scale halos. 

To construct a catalog of cluster-scale halos from the A'- 


body simulation output at z — 0, we used the friends-of- 
friends (FoF) group finder (e.g. Davis et alT) 1985) with a link¬ 
ing length of 0.2 in units of the mean interparticle spacing. 
For each halo we determined the halo center using an itera¬ 
tive technique in which the center of mass of particles within 
a shrinking sphere is computed recursively until a few par- 

[2003, 

by a s 


tides are left inside (e.g. 


Power et al. 


2013|>. Then the halo mass is defined by 


Masaki et al.( 
spherical over¬ 
density method - summing all the particles within a sphere of 
a given overdensity A around the halo center. We constructed 
a catalog that consists of most massive 50 halos from the two 
simulation realizations. Besides the mass threshold, we did 
not employ any other selection criteria such as sphericity or 
the degree of mass distribution complexity. The mean mass 
of the selected halos is similar to the average mass estimated 
from the lensing measurements of 50 Subaru clusters (see Fig- 
ure[2|. Exactly speaking, although the simulated halos are not 
the same in detail as the Subaru clusters, other effects such as 
the intrinsic ellipticities of background galaxies cause much 
larger variations in the lensing profiles as we will show later. 
Hence we believe that the catalog of simulated halos is suit¬ 
able enough for our purpose. 

To test our method as well as to simulate the lensing ob¬ 
servables from the above A'-body simulations, we use the fol¬ 
lowing procedures: 


• 3D mass density profile - We first computed the 

spherically-average mass profile for each simulated 
halo, p(r), where r is the three-dimensional radius from 
the halo center. Then we estimated the NFW parame¬ 
ters, M a and c A for A = 500, by fitting the model NFW 
profile (Eq. ITb to the mass profile, where we weighted 
the simulated mass density profile at a given radial bin 
by the volume of the spherical shell. We stored the best- 
fit parameters c^® 1 ) for each of the 50 halos. 

• 2D lensing profiles - To simulate the lensing profiles 
due to a simulated halo, we use the dark matter (A'- 
body) particles inside or surrounding the halo in the 
simulation output. We estimated the shear profile of 
each halo by projecting the A'-body particles along the 
line-of-sight direction: 

A E(/0 = <Z> «*)-£(*). (19) 


Here we chose the ^-direction of simulation realization 
for the projection, and R is the projected radius from the 
halo center in the xy-plane (the plane perpendicular to 
the projection direction). (X) (< R) is the averaged sur- 
face mass density within a circle of radius R, and £(/?) 
is the averaged surface mass density over the annulus of 
radius R. In this projection calculation, we used a cu¬ 
bic region containing the halo at the center, whose side 
length is 20 h 'Mpc. Since the shear field arises from 
the tidal field around a halo, the constant mass density 
field or the mass density field beyond the cubic region 
causes a negligible contribution to distortion of back¬ 
ground galaxies. We checked that the cubic box is large 
enough for the range of radii we consider. We included 
the reduced shear correction to compute the distortion 
profile of the halo, which is a direct lensing observable: 


AT.(R) = 


AX(fl) 

i-m)/z cr ( Zh z s y 


( 20 ) 


where X cr (zi,z s ) ' s the lensing efficiency. In doing so 
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Mass [10 14 /i _1 Af Q ] 

Figure 2. The distribution of halo masses in 50 clusters, taken from the N- 
body simulations (see text for details) or estimated based on the X-ray ob¬ 
servables (here we used the hydrostatic equilibrium method). 


we assign the source and cluster redshifts of each of the 
50 Subaru clusters to each simulated halo one-by-one 
in descending order of halo masses, where we used the 
hydrostatic equilibrium mass of X-ray observables in 
this matching. The assignment of X cr becomes relevant 
when we will include the effect of background shape 
noise in the Subaru data on the simulated lensing sig¬ 
nals of /V-body halos. We estimated the NFW profile 
parameters, (^ooc’ c 500 r^ by fitting the NFW lensing 
profile (Eq. p~4] > to the above simulated profile. In this 
fitting we weighted the lensing profile at each radial bin 
by the area of radial annulus. We stored the distortion 
profile, the lensing efficiency function, E cr (zi,z s ), and 
the best-fit NFW parameters (M^^ 1 , for each of 

the 50 simulated halos. 


In Figure [2] we compare the distribution of halo masses of 
simulated halos (upper panel) and 50 clusters (lower), where 
we used the best-fit NFW mass for the simulated 2D lensing 
profile for each simulated halo and used the X-ray hydrostatic 
equilibrium mass for each cluster. The figure shows that the 
simulated halos cover the similar range of halo masses as in 
the Subaru clusters. 

Figure [3] compares the best-fit NFW parameters, Msooc 
and C 5000 estimated by fitting the NFW model to the three- 
dimensional mass density profile or the two-dimensional lens¬ 
ing distortion profile for each of the 50 simulated halos. Even 
if we did not include any effect of measurement errors, the 
NFW parameters inferred from the 3D or 2D fitting gener¬ 
ally differ on individual halo basis. For some halos the 2D 
fitting halo mass is larger than the 3D fitting mass, while the 
2D concentration is smaller than the 3D one. These over- or 
under-estimation would be due to the c-M degeneracy in the 
NFW fitting. These biases might cause a source of systematic 
errors in estimating the NFW parameters from the lensing ob¬ 
servables and then testing the ACDM simulation predictions, 
e.g. whether or not the c-M scaling relation inferred from 
the lensing observabl es is consistent with the /V-body simula- 
tion predictions (e.g., Okabe et al.|2010a 2013| Umetsu et al. 



Figure 3. Comparison of the best-fit NFW parameters, halo mass and con¬ 
centration, estimated by fitting the NFW model to the three-dimensional mass 
density profile (“3D NFW-fit”) or the two-dimensional lensing distortion pro¬ 
file (“2D NFW-fit”), respectively, for the 50 massive halos in A/'-body simu¬ 
lations of ACDM model, respectively (see § |2.4| for details). The mean halo 
mass and the range of their halo masses of simulated halos are similar to those 
of 50 Subaru clusters. 


2014|. This is not the main purpose of this paper, but would 
be worth to further study. 

In Figure [4] we performed a hypothetical experiment of the 
stacked lensing analyses with and without NFW scaling, us¬ 
ing the 50 simulated halos. Note that we here ignored shape 
noise contribution for simplicity. First, the upper-left panel 
shows the stacked lensing profile as well as the lensing pro¬ 
files of individual halos, without NFW scaling, i.e. based on 
the standard method. Again note that we used a fixed range 
of radial bins, 0.14 < ^/[/i^Mpc] < 2.8 as we will do for 
the real data. The scatters of individual lensing profiles are 
significant over a range of the radii. Each profile shows var¬ 
ious features due to the aspherical mass distribution, in con¬ 
trast to an analytical, spherical NFW profile. Interestingly, 
however, the figure shows that the average profile after stack¬ 
ing appears to remarkably well match the NFW profile; the 
stacked profile and the best-fit NFW profile are almost indis¬ 
tinguishable, on top of each other. For the sake of comparison 
we plot the amplitudes and the radius relative to the best-fit 
NFW model of the stacked profile as in Figure IT] the best-fit 
parameters are M]?q 0c =* 4 x 10 14 h~ l M Q and c!?q 0( . =* 2.50. 
These numbers are compared to the averages of their underly¬ 
ing true values: (M® c fit ) =* 4.38xl0 14 h~ l M Q , (c®f) =< 2.76 
or (M^f) =* 4.56 x 10 14 r'M G and (c“f) =* 2.57. Thus 
the stacked lensing tends to underestimate the true mass, con¬ 


2005 j ;MenegheftT et al. |20 10] Becker & kravtsov||201 1 


firming the claims in the previous work (Mandelbaum et al. 

Kra 


den Bosch et al.|2013||Meneghetti et al.|2014). 

The other three panels show the results with NFW seal- 
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Figure 4. Simulated lensing profiles for 50 massive halos in the TV-body simulation (see § |2.4| for details). Upper-left panel : The thin curves show the lensing 
profiles for each of the 50 halos, and the bold black curve is the stacked lensing profile without NFW scaling. The blue curve is the best-fit NFW profile to the 
stacked profile. As in Figure^ we plot these profiles in terms of the “scaled” amplitude, A£ + /(2p c r s ) bf-NFW , as a function of the “scaled” projected radius, 
R/r Jf-NFW, where we used the best-fit NFW parameters of the stacked profile (blue curve). Hence both the quantities in the x- and y-axes are dimension-less, 
and can be directly compared to other panels. The lower plot in each panel shows the fractional difference of each profile relative to the best-fit NFW profile. The 
other three panels show the lensing profiles for the same halos when implementing the NFW scaling analysis for each halo or for the stacked analysis (Eq. 0 
Note that the same color curves in the different panels correspond to the same halo. Upper-right panel. The lensing profile when using the NFW parameters 
of three-dimensional mass profile for each halo, (M^ 0 ^ FW_fit , c 5 oo^ FW_fit )’ NFW scaling analysis. Lower-left panel. The results when using the NFW 
parameters of two-dimensional lensing distortion profile, (^ffoo^ IFW_fit 5 c^Q 0 ^ FW_fit ). Lower-right panel. Similar to the upper-right panel, but using the best-fit 
halo mass of each halo and using the halo concentration inferred from the scaling relation, ca = ca(Ma;z) in DK15. In these three panels, the blue curve is not 
a fit, but rather the NFW prediction itself, / NFW (x) (Eq .[ 5 }, where we included the reduced shear correction from the best-fit NFW model to the stacked lensing 
profile in the upper-left panel. Note that, for all the results, we ignored effects of measurement errors such as intrinsic ellipticities of background galaxies. 
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ing implementation. The lensing profiles of individual ha¬ 
los or the stacked lensing profile are estimated by summing 
the “scaled” amplitude of lensing distortion in each of the 
“scaled” radial bin relative to the NFW predictions of each 
halo (see Eq. p~5|). The different panels are the results when 
using the bestTrt NFW parameters of 3D mass density profile 
for each halo, the NFW parameters of 2D lensing profile, or 
the halo mass of 3D profile, but using the concentration pa¬ 
rameter inferred from the scaling relation, csoor = clMsooc) 
in DR14, respective!)]^] The lower-right panel is closest to 
our main results using the Subaru and X-ray data. All the 
three panels clearly show that the NFW scaling significantly 
reduces the scatters of individual lensing profiles relative to 
the NFW prediction, compared to the standard lensing anal¬ 
ysis (the upper-left panel). The scatters appear to be small¬ 
est when using the best-fit NFW parameters of the 2D lens¬ 
ing profile. Comparing the upper-right and lower-right pan¬ 
els clarifies how the scatters are enlarged due to the lack of 
halo concentration knowledge on individual halo or shows the 
degradation by ignoring the intrinsic scatters of halo concen¬ 
tration in different halos. Encouragingly the figure shows that, 
even without knowledge on the concentration of each halo, the 
NFW scaling can reduce the scatters compared to the upper- 
left panel. However it should be noted that the stacked lensing 
profile shows a sizable deviation from the NFW profile (the 
blue solid curve), compared to the upper-left panel. In sum¬ 
mary these results justify our approach of using the scatters of 
individual lensing profiles in order to test the universality of 
mass density profile in 50 clusters. 

3. APPLICATION 

In this section we apply, as a proof of concept, the method 
developed in the preceding section to the Subaru data for a 
sample of 50 massive clusters. 


3.1. The cluster sample 

3.1.1. Subaru weak lensing data 

For the weak lensing measurements, we use the shape cata¬ 
log of galaxies for the 50 clusters, used in the published work 
of Okabe et~aI7| ( 201 3] >. This is the older version of shap e cata¬ 
log, derived as a part of the LoCu SS collaboration (see Okabe 
|et al.|2010a| Martino et al.|2014l for details). In brief, the 50 
cluster samp e was selected fr om the ROSAT All Sky Survey 
catalogs (jEbeling et al.||1998| |2000| |Bohringer et al.||2004 1 
that satisfy the criteria given as L x [0 .1 - 2.4keVJ/£(z) 27 > 
4.2 x 10 44 erg s -1 , 0.15 < z < 0.30, n H < 7 x 10 20 cnT 2 , and 
-25° < 6 < +65°, where E(z) = H(z)/Hq is the normal¬ 
ized Hubble expansion rate. The criteria on the redshift range 
and the declination are adopted in order to have a sufficiently 
high elevation of these clusters from the Subaru telescope and 
to have an entire coverage of the virial region of these clus¬ 
ters with the field of view of the Subaru Suprime-Cam camera 
( Miyazaki et al.|2002j >. 

All the clusters were observed by Subaru, with two pass- 
bands at least: i or I c data, which was used for the weak lens¬ 
ing analysis in [Okabe et al.| (|2013|l, and the bluer-passband 
data, V or g data. For this paper, we take the position of 
brightest cluster galaxy in each cluster as the cluster center. 


Okabe et al. (2010a) carefully studied a possible miscenter- 
ing effect by comparing the lensing signals of various center 
proxies such as the X-ray peak, and concluded that the mis- 
centering, even if exists, should be well within 100 kpc in 
radius (more exactly, within about 50 kpc in our estimate), 
which is inside the minimum radius used in this paper. 

An important systematic effect in the weak lensing mea¬ 
surements is a possible residual uncertainty in estimation of 
source galaxy redshifts, mainly limited by the two passband 
data alone. |Okabe et al.| ( |2013] l developed a method of making 
a secure sample of background galaxies, which is selecting 
galaxies with color sufficiently redder than the red-sequence 
of early-type galaxies in each cluster region. In other words, 
they found that it is very difficult to select “blue” background 
galaxies from the two passband data alone or such blue galax¬ 
ies always appear to be contaminated by foreground or mem¬ 
ber (therefore unlensed) galaxies. However, this selection is 
conservative and leaves only a small number of galaxies in 
the sample so as to ensure less than 1% contamination or di¬ 
lution effect on the lensing signal, even if exists: the mean 
number density of galaxies is about 5 arcmiiH 2 , a factor 4 
or 5 smaller than the number density of all the galaxies for 
which weak lensing analysis is usable in the original i- or I c - 
band catalog. Hence the measurement errors of weak lensing 
signals are dominated by the shape noise, which justifies that 
we ignore the error contribution of projection effects due to 
different structures along the same line-of-sight to the clus¬ 
ter. The mean redshift of background galaxies in each cluster 
was estimated by matching color of the selected background 
galaxies to the COSMOS catalog. Since all the clusters are 
at low redshift z ~ 0.2 and the deep Subaru data typically 
probe galaxies at z ~ 0.8, the lensing efficiency has a weak 
dependence on source redshift and a possible residual uncer¬ 
tainty in the source redshift would not be large and should be 
less than a 10% change in the lensing amplitude even if exists 
(see § 5.7.2 in Okabe et al.||2010a i. We should also keep in 
mind an additional uncertainty due to the sample variance in 
the COSMOS calibration catalog, which refers a possible dif¬ 
ference in the populations of source galaxies in between the 
COSMOS and cluster regions. 

Again note that the purpose of this paper is to give a proof 
of concept of the novel cluster lensing measurement method, 
so the results we will show below is based on the catalog of 
Okabe et al.|(|20131). See Okabe et al. (2015) for the improved 
results of weak lensing measurements based on a more careful 
treatment of shape measurement and photo-z uncertainty. 


3.1.2. X-ray observables: hydrostatic equilibrium mass, gas mass 
and gas temperature 

All the 50 clusters were observed by the X- ray satellites, 
XMM-Newton or/and Chandra (Zhang et akj 12010, Martino 
et al.||2014 i. In this paper, we use the X-ray observables in 
Martino et al.| ( |20i4| to infer the halo mass for each of the 50 
clusters, which was estimated based on either or both of the 
XMM and/or Chandra data. In the following we will use two 
mass estimates: the mass estimated based on the hydrostatic 
equilibrium assumption (hereafter HSE for simplicity) and the 
self-similar scaling relation of gas mass (M gas ) with the halo 
mass: 


6 For the measurement we used the X-ray observables to infer the halo 
mass of each cluster. Here we assumed that the X-ray observables are sensi¬ 
tive to the inner region of each cluster or relatively less sensitive to the projec¬ 
tion effect than in weak lensing. Hence we assume that the X-ray observables 
gives a proxy of the halo mass via the 3D profile. 


Hydrostatic equilibrium (HSE) mass - Martino et al. 


(2014) carefully developed a method of estimating the 
HSE mass of each cluster by combining the surface 
brightness and temperature profiles, measured from the 
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redshift z redshift z 

Figure 5. The 50 clusters used in this paper, in the cluster mass (Msooc ) and redshift plane. The mass estimate of each cluster was taken from Martino et al.| 
12014) (see also Tablefl), derived based on the Chandra and/or XMM X-ray data. The left panel is the mass estimate based on the hydrostatic equilibrium (HSE) 
assumption, while the right panel shows the results derived using the scaling relation of X-ray gas mass with halo mass (Eq. |22) , respectively. The errorbars of 
each cluster are also taken from Martino et al.1)2014) (for M gas , 5 oo we propagated the errors of gas mass). Comparing the two panels reveals that the same cluster 
(symbols at the same x-axis value) generally has different mass estimates and errorbars. 


Chandra and/or XMM data: 


M(< r) = — 


kT g (r)r 

Gpm 


p 


d\np g (r) dlnT g (r) 
d In r dlnr 


( 21 ) 


where T g (r) and p g (r) are the three-dimensional radial 
profiles of gas temperature and density, respectively. 
We will use the halo mass estimate for A = 500 over¬ 


density, M 5 oor, in Table 2 of Martino et al. (2014). For 
some of the clusters, the mass estimates were derived 
for both the Chandra and XMM data. If the two mass 
estimates are available, we use the XMM- derived mass 
because the mass accuracy is better than that of the 
Chandra- based estimate. Note that the two estimates 
are consistent with each other within the errorbars, We 
use the XMM- based mass for 32 clusters, and use the 
Chandra- based mass for the remaining 18 clusters. 

M gas derived mass - The direct X-ray observables are 
the gas mass and temperature. If non-gravitational pro¬ 
cesses are not significant for cluster evolution, the mass, 
temperature, size and other properties of galaxy clus¬ 
ter follow self-similar scaling relations ( |Kaiser|l986| ). 
The ratio of the total matter and gas masses in a clus¬ 
ter region is expected to follow the scaling relation: 
M a (< r) oc M gas A (< r). For the interior g as mass, we 
will use the M gaSi 500 value in Table 3 of Martino et al. 
( j2014[ ) for each cluster. For the normalization factor, 
we here simply employ the cosmic m ean value that is 
inferred from the latest Planck result (Planck Collabo- 
ration et al.|2015) : 


M^qq c QmO 4f gaS) 500 


10 I4 M o Q b0 10 14 M 0 

/ -4f gas> 5O0 


: 11.6 X 


\ 10 14 /r 3 / 2 M F 


( 22 ) 


where we took the best-fit values of Obo h 2 , Q m o h 2 and h 
in Table 3 of [Planck Collaboration et aTTj ( |2015) to com¬ 
pute the normalization constant. The unit /U^ 2 of gas 


mass is from the fact that the gas mass estimate from 
X-ray observables has the /(-dependence. Note that the 
overdensity radius rsooc used for the interior mass def¬ 
inition is from the total mass profile derived from the 
HSE assumption, Eq. (21 1 . In this sense, exactly speak¬ 
ing, this treatment is not self-consistent. Comparing the 
above normalization constant with Figure 2 in Okabe| 


et al. (2014) shows that our model is within a range 
of the normalization constants implied from observa¬ 
tions. However, a precise determination of the normal¬ 
ization constant is not our primary purpose, and the 
above choice can be considered as a working exam¬ 
ple. We will below study how variations in the above 
scaling relation change the weak lensing measurements 
with NFW scaling. 

Table[l]gives a summary of the above X-ray observables: the 
HSE mass and the gas mass for each. Figure[5]shows the dis¬ 
tribution of 50 clusters in the plane of halo mass and redshift. 
The two proxies give a different estimate of halo mass on in¬ 
dividual cluster basis and the errorbars quoted are also differ¬ 
ent. The mean mass of 50 clusters (without lensing weights), 
(Msooc) /[10 14 /r'M o ] = 4.42 or 3.82 for the HSE or gas mass 
proxy, respectively. 

3.2. The stacked lensing analysis of 50 clusters with and 
without NFW scaling 

First of all, in Figure [6] we show the stacked lensing profile 
of 50 clusters, without NFW scaling, for the sake of compari¬ 
son with the following results. This result reproduces Figure 3 
in |Okabe et al.| ([2013]). We employed 32 logarithmically- 
spaced bins over the radial range of 0.14 < R/[h l Mpc] < 
2.8. As given by Eq. ( fl2| ), we estimated the representative 
value of each radial bin by averaging the centric-radii of back¬ 
ground galaxies in the annulus, and therefore the neighboring 
bins are, exactly speaking, not equally spaced, although the 
difference is very small after the average of 50 clusters. The 
cumulative signal-to-noise ratio is significant: S/N ^ 34.5. 
From the fitting to an NFW profile, we find the best-fit param¬ 
eters, M 500 C = (4.0±0.1)xl0 14 /i _I M G and C 500 C = 2.8+0.3, re- 
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Table 1 

X-ray and Lensing Observables of 50 Clusters 


Cluster 

redshift 

''500c- [Mpc] 

X-ray data 

Telescope M gas ^oo [10 14 M o ] 

M 500c [10 14 M q ] 

M hse M Mm 

.500 r .500r 

Lensing observables 
(S/A)wl (rf 2 ) HSE W 2 ) M « as 

A2697 

0.232 

1.20 ±0.04 

XMM 

0.880 ± 0.037 

6.29 ± 0.65 

5.98 ± 0.25 

6.42 

4.47 

4.62 

A68 

0.255 

1.40 ±0.20 

XMM 

0.903 ±0.135 

10.44 ± 0.86 

6.13 ±0.92 

6.24 

13.62 

7.52 

A2813 

0.292 

1.25 ±0.10 

XMM 

1.010 ±0.092 

6.32 ± 0.69 

6.86 ±0.63 

5.45 

4.66 

4.28 

Al 15 

0.197 

0.89 ± 0.07 

XMM 

0.546 ± 0.089 

3.26 ±0.21 

3.71 ±0.60 

3.61 

8.20 

8.46 

A141 

0.230 

1.02 ±0.12 

Chandra 

0.550 ± 0.084 

3.97 ± 1.37 

3.74 ± 0.58 

5.35 

6.01 

6.36 

ZwC10104 

0.254 

0.76 ± 0.01 

Chandra 

0.235 ± 0.004 

1.67 ±0.07 

1.60 ±0.03 

3.55 

5.46 

5.50 

A209 

0.206 

1.15 ±0.07 

XMM 

0.972 ± 0.094 

5.45 ±0.18 

6.60 + 0.64 

9.14 

21.84 

17.37 

A267 

0.230 

1.17 ± 0.12 

Chandra 

0.703 ± 0.094 

5.97 ± 1.84 

4.78 ± 0.64 

5.37 

7.01 

7.47 

A291 

0.196 

0.94 ± 0.06 

XMM 

0.391 ±0.031 

2.92 ± 0.56 

2.66 ±0.21 

4.92 

8.61 

8.92 

A383 

0.188 

1.01 ±0.08 

XMM 

0.425 ± 0.036 

3.25 ± 0.45 

2.89 ± 0.24 

5.74 

9.11 

9.80 

A521 

0.248 

1.25 ±0.04 

XMM 

1.230 ±0.050 

7.05 ± 0.34 

8.36 ± 0.34 

5.98 

4.66 

5.23 

A586 

0.171 

1.09 ±0.08 

Chandra 

0.600 ± 0.049 

4.42 ± 0.90 

4.08 ± 0.33 

5.58 

10.90 

11.47 

A611 

0.288 

1.20 ±0.06 

Chandra 

0.612 ±0.039 

6.80 ± 1.08 

4.16 ±0.26 

5.98 

8.24 

12.30 

A697 

0.282 

1.50 ±0.10 

Chandra 

1.490 ±0.127 

13.14 ±2.74 

10.12 ±0.86 

4.89 

20.89 

12.84 

ZwC10857 

0.235 

0.85 ± 0.04 

XMM 

0.351 ±0.017 

2.33 ±0.23 

2.38 ±0.12 

2.45 

4.53 

4.60 

A750 

0.163 

0.97 ±0.12 

Chandra 

0.406 ± 0.061 

3.17 ± 1.22 

2.76 ±0.41 

6.82 

16.54 

18.17 

A773 

0.217 

1.21 ±0.11 

XMM 

0.907 ±0.102 

7.77 ± 0.94 

6.16 ±0.69 

9.03 

7.40 

11.40 

A781 

0.298 

1.13 ±0.17 

XMM 

0.783 ±0.142 

6.72 ± 0.45 

5.32 ±0.96 

4.27 

28.93 

25.32 

ZwC10949 

0.214 

0.93 ± 0.02 

Chandra 

0.315 ±0.010 

2.90 ± 0.20 

2.14 ±0.07 

5.50 

15.34 

16.79 

A901 

0.163 

0.79 ± 0.06 

XMM 

0.208 ± 0.020 

1.68 ±0.39 

1.41 ±0.14 

5.58 

11.98 

13.43 

A907 

0.167 

1.08 ±0.06 

XMM 

0.623 ± 0.042 

5.17 ±0.64 

4.23 ± 0.29 

7.71 

11.32 

14.53 

A963 

0.205 

1.14 ±0.07 

XMM 

0.616 ±0.050 

5.60 ±0.71 

4.18 ±0.34 

7.49 

13.40 

17.46 

ZwC11021 

0.291 

1.26 ±0.05 

XMM 

1.080 ±0.048 

6.82 ±0.14 

7.34 ±0.33 

6.78 

8.76 

9.23 

A1423 

0.213 

1.18 ± 0.10 

Chandra 

0.711 ±0.095 

6.02 ± 1.53 

4.83 ± 0.65 

4.58 

9.89 

6.95 

A1451 

0.199 

1.36 ±0.11 

XMM 

1.050 ±0.102 

8.97 ±2.18 

7.13 ±0.69 

8.25 

6.30 

5.44 

RXCJ1212 

0.269 

0.76 ± 0.05 

XMM 

0.196 ±0.012 

1.67 ±0.31 

1.33 ±0.08 

3.00 

9.34 

8.86 

ZwC11231 

0.229 

1.23 ±0.08 

Chandra 

0.828 ± 0.078 

6.82 ± 1.29 

5.63 ±0.53 

5.09 

17.35 

16.05 

A1682 

0.226 

1.24 ±0.18 

Chandra 

0.764 ±0.137 

7.35 ± 3.06 

5.19 ±0.93 

7.30 

6.12 

8.03 

A1689 

0.183 

1.52 ±0.07 

XMM 

1.290 ±0.059 

11.98 ± 1.94 

8.76 ± 0.40 

9.42 

11.25 

16.82 

A1758N 

0.280 

1.38 ±0.07 

Chandra 

1.220 ±0.062 

10.21 ± 1.54 

8.29 ± 0.42 

3.10 

13.65 

10.03 

A1763 

0.228 

1.33 ±0.11 

XMM 

1.230 ±0.135 

6.60 ± 0.56 

8.36 ±0.92 

6.75 

9.58 

7.09 

A1835 

0.253 

1.57 ±0.11 

XMM 

1.550 ±0.120 

14.04 ± 1.27 

10.53 ± 0.82 

5.96 

14.77 

13.27 

A1914 

0.171 

1.38 ±0.08 

XMM 

1.160 ±0.073 

8.08 ± 1.00 

7.88 ± 0.50 

4.83 

9.22 

8.88 

ZwC11454 

0.258 

1.06 ±0.10 

XMM 

0.578 ± 0.060 

3.65 ± 0.42 

3.93 ±0.41 

3.18 

3.94 

4.13 

A2009 

0.153 

1.29 ±0.13 

Chandra 

0.708 ± 0.082 

7.33 ± 2.47 

4.81 ±0.56 

4.91 

5.79 

5.14 

ZwC11459 

0.290 

1.08 ±0.25 

XMM 

0.675 ±0.187 

5.65 ±0.36 

4.59 ± 1.27 

3.31 

4.16 

3.73 

RXCJ1504 

0.215 

1.47 ±0.35 

XMM 

1.300 ±0.337 

10.93 ± 0.82 

8.83 ± 2.29 

4.55 

7.70 

6.23 

A2111 

0.229 

1.17 ± 0.14 

Chandra 

0.719 ± 0.110 

5.99 ± 1.89 

4.88 ± 0.75 

4.83 

13.10 

12.57 

A2204 

0.152 

1.49 ±0.08 

XMM 

1.280 ±0.082 

10.66 ± 1.72 

8.70 ±0.56 

6.48 

11.20 

9.16 

A2219 

0.228 

1.75 ±0.11 

XMM 

1.882 ±0.216 

14.35 ± 2.04 

12.79 ± 1.47 

7.51 

6.21 

6.06 

RXCJ1720 

0.164 

1.23 ±0.11 

XMM 

0.771 ±0.083 

6.97 ± 0.68 

5.24 ±0.56 

3.80 

11.47 

7.02 

A2261 

0.224 

1.22 ±0.12 

Chandra 

1.000 ±0.127 

6.75 ± 1.89 

6.79 ± 0.86 

8.88 

13.44 

13.33 

RXCJ2102 

0.188 

1.00 ±0.06 

XMM 

0.450 ± 0.033 

3.52 ±0.61 

3.06 ± 0.22 

4.04 

13.87 

12.78 

RXJ2129 

0.235 

1.08 ±0.04 

XMM 

0.749 ± 0.037 

4.22 ±0.16 

5.09 ± 0.25 

3.17 

5.10 

6.17 

A2390 

0.233 

1.60 ±0.11 

XMM 

1.700 ±0.088 

13.67 ±2.09 

11.55 ±0.60 

6.30 

8.47 

7.41 

A2485 

0.247 

1.11 ±0.15 

Chandra 

0.558 ± 0.087 

5.32 ±2.08 

3.79 ± 0.59 

4.71 

1.89 

0.52 

A2537 

0.297 

1.19 ± 0.10 

XMM 

0.739 ±0.081 

7.20 ±0.73 

5.02 ±0.55 

5.05 

8.73 

9.81 

A2552 

0.300 

1.25 ±0.09 

Chandra 

1.020 ±0.094 

7.81 ± 1.64 

6.93 ± 0.64 

3.98 

11.30 

10.85 

A2631 

0.278 

1.20 ±0.09 

XMM 

1.030 ±0.088 

8.51 ±0.98 

6.93 ± 0.64 

4.61 

18.05 

16.23 

A2645 

0.251 

1.15 ± 0.18 

Chandra 

0.541 ±0.117 

5.98 ± 2.59 

3.68 ± 0.79 

6.59 

19.54 

19.62 


Note. 
A = 500. 


— The X-ray observables (75000 an ^ Af gas, 500 ) taken from Tables 2 and 3 of 
^ 500 c * s mass estimate based on the hydrostatic equilibrium, and M^ 


Martino et al. 


2014 : rsooc is the radius for the interior overdensity 

Meas is the total 


gas,500 is the gas mass interior to 7*500 (see text for details). M 5QOc 

mass interior to rsooc assuming the simple self-similar scaling relation given by Eq. \22) . The mean mass of 50 clusters (M 5 00c) /[10 14 /z _1 M o ] = 4.4 2 or 3.82 for 
t he HS E and gas mass cases, respectively. The last three columns are the lensing observables that are computed from the lensing measurement of|Okabe et al.| 
ra. (S/A0 wl is the total signal-to-noise ratio of lensing distortion measurement for each cluster over the 8 radial bins in the r ange 0.14 < R/\h *Mpc] < 2.8. 
^2)HSE or ^2^M gas j s the deviation of the lensing distortion profile compared to the NFW prediction, defined by Eq. 0 or jl8[ for each cluster. 
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Figure 6. Upper panel : The stacked distortion profile measured from 50 
Subaru clusters, based on the standard stacked lensing analysis (Eq. [9j, i.e. 
without NFW scaling implementation. We employed the 32 logarithmically- 
spaced bins over a fixed range of radii, 0.14 < 7?/[/i _1 Mpc]< 2.8, for all the 
50 clusters. The errorbar at each bin is computed from Eq. 0 assuming that 
the statistical noise is dominated by the intrinsic ellipticities of background 
galaxies. The solid curve is the best-fit NFW model, which is specified by the 
best-fit parameters Msoo c = (4.0 + 0.1) x 10 l4 h~ l M Q and C 500 C = 2.8 ± 0.3. 
The bestTit model reproduces the results in jOkabe et al.| {2013) within the 
errorbar^] For the sake of comparison with the following figures, we plot 
the distortion profile relative to the best-fit NFW model, as a function of the 
radius relative to the scale radius of the best-fit NFW model. Note that the 
representative value of each radial bin is estimated from t he av erage of radii 
of background galaxies that reside in the annulus (see Eq. |12) . The reduced 
chi-square for the best-fit model is^ 2 /d.o.f = 22.5/(32 - 2). Lower panel. 
Similar to the above panel, but for the 45°-rotated components of background 
galaxy ellipticities. 

JOkabe et al.|(|2013) used the physical coordinates in the stacking analysis, 
and this leads to a slight change in the best-fit NFW parameters. 

spectively. The reduced chi-square is^f 2 /d.o.f = 22.5/(32-2). 
Thus the results show that, even if the X-ray inferred masses 
differ from each other by up to a factor of 10, the stacked pro¬ 
file is so remarkably well fitted by the NFW model. This ap¬ 
pears to be consistent with what we found from the test using 
the simulated halos in Figure [4] 

We now move to the main results of this paper. We em¬ 
ploy the following procedures to implement the NFW scaling 
analysis of weak lensing measurements: 

(1) NFW scaling of galaxy ellipticities and radial bins - 
First, we employ, for the fl-th cluster (a = 1,2, • • ■ , 50), 
the halo mass inferred from the X-ray observables, ei¬ 
ther HSE or gas mass (see § |3.1.2[ ). We then use 
the c-M relation in DK15 to infer the halo concen¬ 
tration for the cluster. Using the X-ray inferred pa¬ 
rameters, Mf 00(a) and c? (1(|(a) , we compute the expected 
lensing amplitude and the NFW scale-radius, 2[p r r v ] (a| 
and r s ( a ), respectively, in order to “scale” the ampli¬ 
tude of galaxy ellipticities as well as the radius for 
the Sfl-th background galaxy in the a-th cluster region: 
e+(v„i * <2+(.v„j/ [2p,rf( al and R( a )s a > x = R( a )s a / r s{a)- 
Thus, even if we use the same background galaxies 
over a fixed range of radii, 0.14 < /?/[/i _1 Mpc] < 2.8, 
this NFW scaling leads to different ranges of the scaled 



Figure 7. The expected, differential signal-to-noise ratio, d(5 /N), at each of 
the 8 logarithmically-spaced radial bins in the range 0.14 < /?/[/i -1 Mpc] < 
2.8, for each of 50 Subaru clusters. We computed the d{S/N) value as fol¬ 
lows. For the expected signal, we used an analytical NFW profile for each 
cluster assuming the X-ray HSE mass and the halo concentration inferred 
from the c-M relation in DK15. To compute the statistical noise in each bin, 
we used the real Subaru data of background galaxies (their distribution on the 
sky, the intrinsic shapes and the lensing weights) in each cluster region. Most 
of the data points, about 79% of 400 data points (400 = 50 x 8), are expected 
to have the d(.S /N) values greater than unity. The same color symbols at dif¬ 
ferent radial bins correspond to the same cluster. No te that the representative 
value of each radial bin is computed from Eq. 0 taking into account the 
radii and weights of background galaxies, which causes variations in the rep¬ 
resentative values especially for the small radii, even if we work on the fixed 
range of 0.14 < R/[/t _1 Mpc] < 2.8. 

radii, x, for different clusters. The different amount of 
radial scaling requires a careful treatment of the radial 
binning, especially when comparing the lensing distor¬ 
tion profiles with and without NFW scaling. In the fol¬ 
lowing we use the different binning schemes depending 
on either case studying the stacked lensing profile or 
studying the scatters of 50 cluster lensing profiles rela¬ 
tive to the NFW prediction, which are summarized by 
the procedures (2a) and (2b) below. 

(2a) Stacked lensing analysis with NFW scaling implemen¬ 
tation - As in Figure [61 we will study the stacked dis¬ 
tortion profile of 50 clusters after the NFW scaling of 
each cluster. Similarly to Figure [6] we will use the 32 
logarithmically-spaced bins in the “scaled” radius, x, 
where we used exactly the same background galaxies 
behind the 50 clusters. After stacking 50 clusters, we 
can expect a significant detection of the lensing signal 
at each radial bin, as implied from Figure [6] However, 
the above NFW scaling transforms the original radial 
range to a different range of the scaled radius x for each 
cluster. Hence, the sample of background galaxies in 
each bin of R or x radii differ from each other. Never¬ 
theless, since the stacked lensing has a sufficiently high 
S /N at each bin, we checked that the NFW scaling al¬ 
most conserves the total S /N value (exactly speaking, 
it causes only about 0.5% fractional change). We es¬ 
timate the representative value of each radial bin in a 
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Figure 8. Similar to Figure|6] but the stacked distortion profile of 50 clusters when implementing the NFW scaling analysis: we summed the “scaled” amplitudes 
of background galaxy ellipticities in each bin of the “scaled” radii according to the NFW parameters, halo mass and concentration, inferred for each cluster based 
on its X-ray observables. In the left or right panels, we employed the X-ray inferred mass of each cluster from the hydrostatic equilibrium assumption (HSE) 
or the gas mass, respectively, and then used the halo concentration inferred from the scaling relation c = c(M;z ) in DK15. Note t hat w e used exactly the same 
background galaxies as those for the analysis without NFW scaling in Figure [6] The errorbar at each bin is computed based on Eq. i |16| >. The soli d cu rve in each 
panel is not a fit, but the NFW prediction (y NFW given by Eq.[5j including a small correction due to reduced shear at the small radii (see below Eq. |l5| for details). 
The reduced chi-square is^ 2 /d.o.f = 31.3/32 or 30.7/32 for the HSE or gas mass case, respectively. 
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RXCJ2102 

ZwC10949 

RXJ2129 

A0901 

A2390 

A0907 

A2485 

A0963 

A2537 

ZwC11021 
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x = R/r s 


Figure 9. The difference between the lensing distortion profiles of 50 clusters and the best-fit NFW profile (A£ best fit (/?)) or the normalized NFW profile 
(/NFW(x)) for the weak lensing analysis with or without NFW scaling implementation in the left or right panel, respectively, as in the lower panels of Figure [ 4 ] 
Here we show in the right panel the result when using the gas mass to estimate the halo mas s of ea ch c luster, as in the right panel of Figure [8] To make a fair 
comparison, we show the relative difference to the statistical error at each radial bin (see Eqs.|17|and |18| . Since the lensing profile is noisy on individual cluster 
basis due to the fewer number of background galaxies, we employed the 8 logarithmically-spaced bins in the fixed range of 0.14 < 7?/[/z -1 Mpc] < 2.8 for all the 
clusters as in Figure]?] In addition, we used the same background galaxies in each radial bin before and after the NFW scaling tra nsfor mation, x = R/r s , for each 
cluster so that the differences become identical if we set the model NFW profile AE bf-NFW = / NFW = 0 (see the procedure 2b in § |3.2| for details). Also note that, 
due to the NFW scaling, the fixed radial range in the left panel is transformed to the different range of the scaled radius for different clusters. The same-color 
curves in the two panels correspond to the same cluster, and the bold curve shows, as an example, the result for A781, which has the largest deviation from the 
NFW profile. Since the sum of squares of all the curves gives an estimate to quantify the scatters of 50 cluster lensing profiles relative to the NFW model - we 
call the d 2 value. The NFW scaling yields d 2 = 527.1 or 504.6 for the HSE and gas mass cases, respectively, compared to d 2 = 543.2 for the case without NFW 
scaling (Figure jbj. This corresponds to the improvement A d 2 = d 2 - ^-scaling = (4-0) 2 or (6.2) 2 , respectively. 
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Figure 10. Similar to the previous figure, but the scatters for 25 clus ters, 
which are a half of 50 clusters that have the larger value of A d 2 (Eq. \23\ 
than other 25 clusters on individual cluster basis. Comparing the left ana 
right panels more clearly shows that the NFW scaling reduces scatters of the 
lensing profiles. The cumulative value of A d 2 for 25 clusters, A d 2 — (9.3) 2 . 


similar manner to Eq. m 

(2b) Studying the scatters of l ensin g profiles for 50 clusters 
- As we discussed in § |2.3| we monitor the scatters 
of 50 cluster lensing profiles relative to the NFW pre¬ 
diction in order to address the existence of the univer¬ 
sal NFW profile. To quantify the scatters, we compute 
the df value for either case with or without NFW scal¬ 


ing (see Eqs. 17 and 18 for the definition). In doing 
this, we need to probe tKe “shape” of lensing profile for 
each cluster, and in other words each radial bin needs 
to be in the signal dominated regime on an individual 
cluster basis. Hence, if we take the 32 bins as in the 
stacked lensing analysis, each radial bin suffers from 
the shape noise contamination due to too low number 
density of background galaxies in each bin. To tackle 
this obstacle, we employ 8 logarithmically-spaced bins 
in the range 0.14 < /?/[/i _1 Mpc] < 2.8 for each clus¬ 
ter. Figure [7] shows the expected S /N at each radial 
bin for the 50 clusters. The figure shows that 319 data 
points among 400 points, corresponding to 79% of 400 
data points, are expected to have the S /N value greater 
than unity. Hence the 8 bins seem suitable for our pur¬ 
pose. TablefTlgives the total S/N of each cluster when 
employing the 8 bins. However, the expected lensing 
signal at each radial bin would be still noisy. To avoid 
any artifact arising from the noise dominated bins, we 
transform each of the original bins in R to the corre¬ 
sponding bin in the scaled radius x after the NFW scal¬ 
ing, rather than redefining the radial bins for a fixed 
range of x. With this binning, each radial bin before 
and after the NFW transformation preserves the same 
background galaxies. Hence, this binning method pre¬ 
serves the S /N value in each radial bin as well as the 
total S /N value for each cluster, before and after the 
NFW scaling, as can be found from Eqs. ( [T7| and (fl8 i 
mathematically. As a result, the different clusters cover 
different ranges of the scaled radius x. 



Figure 11. A test of the performance of the NFW scaling analysis in Fig¬ 
ure [8] We randomly assigned the X-ray inferred mass to each clu ster, redid 
the scaling analysis, and then computed the d 2 difference (Eq.[23). The his¬ 
togram shows the distribution of 200 random realizations, which can be com¬ 
pared to our main result shown by the vertical line for either the HSE mass or 
the gas mass case in the upper or lower panels, respectively. All the random 
realizations give a negative value of A d 2 , and any of those does not reproduce 
the measurement value. Compared to the mean and variance of the random 
realizations, the measured A d 2 value is away from the mean at 3.6 and 3.7 cr 
for the HSE and gas mass cases, respectively. 


model to the stacked distortion profile in Figure [6] The re¬ 
duced shear correction is not large (by up to 20% in the am¬ 
plitude at the inner bins) over the range of radii, as can be 
found from Figure|T| The figure shows that the stacked profile 
is in excellent agreement with the NFW prediction, to within 
the errorbars. This agreement supports the existence of NFW 
profile in the clusters, and implies that the X-ray inferred mass 
indeed gives a proxy of the genuine mass for each cluster. To 
be more precise, the reduced chi-square is^f 2 /d.o.f = 31.3/32 
or 30.7/32 for the HSE or gas mass case, respectively, com¬ 
pared to^- 2 /d.o.f = 22.5/(32 - 2) in Figure[6] 

Now we compare the scatters of 50 cluster lensing profiles 
with and without implementation the NFW scaling, quantified 
by the d 2 value (Eqs.|l7|and|l8|, in order to address the exis¬ 
tence of universal NFW profile. By using the above method 
(2b), we find the difference between the df values with and 
without NFW scaling as 

Ad = d~ — <f w _ sca jj n g 

_ ( 543.2 - 527.1 ^ (4.0) 2 , (HSE) 

\543.2 - 504.6 - (6.2) 2 , (M gas -M 500f ) ’ 


Figure [8] shows the stacked lensing profiles after imple¬ 
menting the NFW scaling (the above case 2a), using the halo 
mass proxies based on the HSE assumption or the gas mass, 
respectively. We again note that, to have a fair comparison 
with Figure [6] we have used exactly the same background 
galaxies. The solid curve in each panel is not a fit, but rather 
is the NFW prediction (Eq. [5]), including the reduced shear 
correction 1/[1 - k nfw (x)], where we used the best-fit NFW 


Thus the NFW scaling for both the HSE and gas mass cases 
leads to the smaller t/ 2 -values, meaning the smaller scatters 
of lensing profiles relative to the NFW profile than the scat¬ 
ters without NFW scaling. The smaller d 2 value for the gas 
mass implies that the gas mass gives a better proxy of the 
underlying true masses of clusters (at least the relative mass 
differences between different clusters). Thus the NFW scal¬ 
ing gives about 4- or 6-cr improvement for the HSE or gas 
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Figure 12. As in the previous figure, but show effects of the statistical errors 
in X-ray inferred halo mass of each cluster on the Ad 2 value. Here, we added 
a random scatter to each halo mass by an amount of the quoted errorbar in 
Table ^assuming the Gaussian distribution, redid the weak lensing analysis 
with NFW scaling, and then computed the Ad 2 value for each realization. 



Figure 13. The figure shows the change in Ad 2 (the vertical axis) when using 
the different range of radii (the horizontal axis) used in the analysis. The bold 
lines are the results for our fiducial choice: 0.14 < R/[h ] Mpc] < 2.8 and 
8 logarithmically-spaced bins for all the 50 clusters. The other lines are the 
results when excluding the innermost or outermost bin from the analysis, 
respectively, for the the HSE or gas mass proxy cases, respectively. 


mass, respectively, assuming that the d 2 distribution obeys 
a ^-distribution. We checked that, even if we use 16 bins 
instead of our fiducial 8 bins, the d 2 values themselves get 
enlarged because each bin is more in the shape noise domi¬ 
nated regime, but the d 2 difference, the Ad 2 value, is almost 


unchanged. 

The smaller d 2 value due to the NFW scaling arises from 
two parts: the scaling of lensing profile amplitude (or back¬ 
ground galaxy ellipticities) and the scaling of cluster-centric 
radius. The two scalings are specified by halo mass and 
concentration of each cluster: e+/[2p c r s ] oc M^ 2 c 5 q 0c and 

x = R/r s oc M 50 q c C5ooc, respectively. If we include only the 
scaling of background galaxy ellipticities, without the radial 
scaling, we found <^_ scaling = 512.9 or 526.5 for the HSE or 

gas mass case, which equivalently correspond to Ad 2 - (5.5) 2 
or (4.1) 2 , respectively. That is, the HSE case shows an even 
greater improvement in Ad 2 compared to Eq. (231. On the 
other hand, if we include only the scaling of radius, but with¬ 
out the scaling of galaxy ellipticities, d 2 = 547.1 or 522.3, 
which correspond to Ad 2 - -3.9 or (4.6) 2 , respectively. Thus, 
for the HSE case, the radial scaling does not appear to be ade¬ 
quate, and rather gives a positive Ad 2 . For the gas mass case, 
both the two scalings about equally contribute to the reduce 
in the d 2 value or the significance of the smaller scatters. 

Figure |9] shows the contribution of each cluster to the up¬ 
value, which shows the argument of Eqs. ( [T7] > or (fl8]l at each 
radial bin for each of 50 clusters. The total d 2 value is ob¬ 
tained by summing the square of each curve over the 8 ra¬ 
dial bins and 50 clusters. Table [l] gives the total f/ 2 -value 
for each cluster. The figure shows that, although it looks 
noisy, the NFW scaling reduces the scatters. One might no¬ 
tice some outlier clusters: the clusters, which have top three 
largest ^ 2 _ scalin „ values (see Table IT 1 , are A781, A209, and 
A697 for™ the HSE case, while A781, A2645 and A750 for 
the gas mass case, respectively. For example, the bold curve 
denotes the A781 cluster, which is a more than 3cr outlier 
than the mean. The mass distribution of A781 displays three 
prominent peaks, indicating that the cluster is in the phase of 
ongoing merger (Wittman et al. 2014). 

For further clarification, - Figure flOl shows the scatters for 
25 clusters which are a half of 50 clusters that have the 
larger value of Ad 2 (Eq. 23), i.e. show the better improve¬ 
ment of NFW scaling analysis, than other 25 clusters (see 
Table |T]). The reduce in the scatters due to the NFW scal¬ 
ing is more evident. In this case, the cumulative value of 
Ad 2 for 25 clusters, (A d 2 ) 1/2 = (86.2) I/2 - 9.3. Similarly, 
(Ad 2 ) 1 ' 2 =* (80.6) I/2 =* 9.0 for the HSE mass. The other 25 
clusters yield smaller or even positive Ad 2 values, probably 
due to the complex mass distribution or the inaccuracy in X- 
ray inferred halo masses. 

To draw a more robust conclusion, we make several tests of 
our results. In Figure [TT| we studied how the scatters of 50 
lensing profiles are enlarged if we implement the NFW scal¬ 
ing analysis by randomly assign the X-ray inferred halo mass 
to each cluster (without repeated use of X-ray mass). All the 
200 random realizations have a negative value of Ad 2 , and any 
of the random realizations cannot reproduce a similar positive 
value to the measured A dr (the vertical line) for both the HSE 
and gas mass cases. To be more quantitative, the measured 
value Ad 2 is away at 3.6 and 3.7 cr for the two cases, respec¬ 
tively, compared to the mean and variance of the random real¬ 
ization distribution. These results give another support on the 
existence of NFW profile in the 50 clusters. 

One important source of uncertainties in the method is a 
residual uncertainty in the X-ray inferred halo mass or a pos¬ 
sible effect of intrinsic scatter in the mass scaling relation of 
X-ray observable. Figure [12] shows how the statistical errors 
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of X-ray inferred mass affect the Ad 2 value. To be more pre¬ 
cise, we added a random scatter to halo mass of each clus¬ 
ter assuming the Gaussian distribution with variance given 
by the quoted errorbar of each X-ray mass in Table [T| i.e. 
M'oo (l() = M^ 00(a) + < 5M (a ), treated the shifted mass as its true 
mass, and then redid the NFW scaling analysis. The figure 
shows that adding the random scatter to each cluster tends to 
decrease Ad 2 , implying that the central value of the X-ray in¬ 
ferred mass is indeed closer to the underlying true mass. The 
distribution of Ad 2 is wider for the HSE mass, but this would 
be ascribed to the larger errors of HSE mass than those for 
the gas mass as can be found from Figure [5] Again encour¬ 
agingly, even if adding the random errors tothe gas mass, the 
resulting Ad 2 values are positive, supporting that the gas mass 
is a better proxy of the genuine cluster mass as in Figure [8] 
Although we have used the fixed range of the original co¬ 
moving radius, 0.14 < R/[h 'Mpc] < 2.8, for all the clus¬ 
ters as our fiducial choice, Figure [13] shows how the results 
are changed if excluding the inner- or outermost radial bin of 
8 logarithmically-spaced bins from the analysis. The figure 
shows that, for the X-ray gas mass proxy, excluding the outer- 
or innermost bin degrades the NFW scaling or reduces the Ad 2 
values, suggesting that the wider range of radii is important to 
capture the curvature of the mass profile. On the other hand, 
for the HSE mass case, excluding the outermost bin increases 
the Ad 2 , again implying that the HSE mass estimate might not 
be as accurate to infer the genuine mass as the gas mass and 
involve residual systematic errors. 


3.3. Discussion and Implications 
3.3.1. Comparison with N-body simulations 


Is the NFW scaling results in Figures [9] and Eq. (23 i as 
expected from a viewpoint of ACDM structure formation 
model? To address this question. Figure pA] compares the 
measurement results and the /V-body simulated halos, as in 
Figure P] in a two-dimensional space of the d 2 values with 
and without NFW scaling. To make a fair comparison, we 
included the effect of intrinsic galaxy ellipticities on the sim¬ 
ulation results. To be more precise, (1) we first populated, 
into each region of simulated halos, the background galaxies 
taken from the corresponding Subaru cluster data (matched in 
descending order of halo masses), (2) made a random rotation 
of orientation of each galaxy ellipticity, which erases the co¬ 
herent lensing signal of each Subaru cluster, (3) simulated the 
“observed” galaxy ellipticity of each background galaxy by 
adding both the lensing distortion of simulated halo and the 
intrinsic shape, (4) performed the hypothetical lensing mea¬ 
surements with and without NFW scaling, and (5) computed 
their d 2 values, respectively. To account for the statistical vari¬ 
ance of intrinsic ellipticities, we generated 40 realizations of 
the /V-body simulation results: we redid the d 2 calculations 
after random rotation of background galaxies. For the sim¬ 
ulation results, we consider the three cases similarly to Fig¬ 
ure [4] the lensing analysis with NFW scaling when using the 
best-fit NFW parameters of 3D mass profile, the NFW pa¬ 
rameters of 2D distortion profile, or the best-fit halo mass of 
3D profile, but using the halo concentration inferred from the 
c-M relation, respectively. The third case is closest to what 
we did for the actual data. First of all, the simulation re¬ 
sults without NFW scaling, denoted by the d 2 values in the 
horizontal axis, fairly well reproduce the measurements on 
average, reflecting that the statistical errors in the d 2 value 
are dominated by the shape noise. Also note that the hori- 



d 2 


Figure 14. Comp ariso n of th e measurements and the simulation results for 
the devalues (Eqs. |l7| and |l8| for their definitions). The bigger-size circle and 
square symbols are the measurement results for the HSE and gas mass cases, 
respectively. The other symbols are the simulation results as in Figure [4] but 
including the effect of intrinsic ellipticities of background galaxies that are 
taken from actual Subaru data of each cluster region (see text for details). The 
tilted triangle and hexagon symbols are the simulation results when using the 
best-fit NFW parameters of 3D mass profile for each halo, the NFW param¬ 
eters of 2D projected lensing profile or the best-fit mass of 3D mass profile, 
but the concentration inferred from the scaling relation = ca(Ma;z), re¬ 
spectively (see Figure [4]for details). The third case is intended to mimic what 
we did for the actual measurements. We show the simulation results for 40 
realizations of background galaxy ellipticities. Note that, for each realization, 
we computed the three simulation results; each of the triangle and hexagon 
symbols has the same d 2 value in the horizontal axis, but different d 2 .. 

J ’ w—scaling 

values in the vertical axis. For comparison, the orange-color star symbol de¬ 
notes one particular realization that has a similar d 2 value to the measurement 
for no NFW scaling case (the horizontal axis). The two star symbols in the 
left-lower comer are the results when using the same realization of back¬ 
ground galaxies as in the orange-color star symbol, but using the analytical 
NFW profiles for the d 2 calculations. Note that the arrow in the lower-left 
comer denotes the simulation result that is below the plotted range. 


zontal spread of the simulation realizations is roughly given 
by vd? - V550 ^ 23. However, all the simulation results 
with NFW scaling, d 2 , in the vertical axis, are system- 

atically smaller than the measured values. Thus this disagree¬ 
ment suggests that we do not properly consider some effects 
inherent in the measurements on the simulation results. For 
comparison, the star symbols show the results when using an¬ 
alytical NFW halos to compute the d 2 values where we used 
the X-ray proxy masses for the HSE or the gas mass to com¬ 
pute the NFW lensing profile of each cluster. The difference 
between the analytical NFW halos and the simulation results 
is due to the complexity of mass distribution in the simulated 
halos, such as asphericity, substructures, and the scatters of 
halo concentration. 

A possible source to reconcile the difference between the 
measurements and the simulation results in Figure 14 is an 


additi onal error or in trinsic sc atter in the X-ray inferred halo 
mass (|Stanek et al.||2010| |Okabe et al.||2010b|). Figure [T5] 
addresses this question. The left panel shows how adding 
a scatter to each mass of simulated halos, parametrized by 
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Figure 15. Sim ilar to the previous figure, but here included the effects of mass scatter of each halo on the simulation results. The hexagon symbols are the same 
as in Figure p~4] For each realization of background galaxies, we added a random mass scatter to each halo, simulated the lensing analysis with NFW scaling by 
treating the shifted mass as the true mass, and then computed the <^_ scalin value. Adding the halo mass scatters tends to degrade the NFW scaling results or 
preferentially causes an up-scatter of each simulation result in this two-dimensional space. Left panel'. The simulation results when adding the Gaussian mass 
scatters by the fractional errors of o"i„m = 0.1,0.2 or 0.3, respectively. The arrows in the upper horizontal axis denote the case that the simulation results are 
outside the range shown in this plot. Right panel: Similar to the left panel, but the results when adding a random mass scatter to each simulated cluster assuming 
the fractional error proportional to the quoted error of the gas mass proxy relation (assigned to the simulated halo); cn n M (0) = er^ //V/^ or a factor 2 or 3 bigger 
one. 


the fractional variance cr lnM = cr(M)/M = 0.1,0.2 or 0.3, 
degrades the ^_ scaling values for the NFW scaling analysis. 
More precisely, we randomly generated a mass scatter 6M 
for each halo assuming the Gaussian distribution with vari¬ 
ance cr\ n M, added the scatter to each halo mass as given by 
M\ a) = -^500(0’) + and then computed the ^J,_ scalino value 
by treating the shifted mass M' as the true mass of each sim¬ 
ulated halo. For the sake of comparison, we used the same 
40 realizations of background galaxies as in Figure [14] and 
therefore the degradation is solely due to the mass scatters. 
Note that, for each realization of background galaxies, adding 
the halo mass scatters changes only the <^_ scaling value in the 
vertical axis. The figure shows that the halo mass scatters 
generally degrades the NFW scaling result or equivalently en¬ 
large the di sca | in „ value. However, only the additional errors 
of cr lnM ~ 0.2-0.3 can reproduce the measurement result for 
the gas mass proxy. This might imply that the X-ray halo mass 
involves an unknown, systematic error or intrinsic scatter. 

As an alternative test, the right panel of Figure [T5] shows 
the effects of the quoted errorbars in the X-ray inferred halo 
masses. Here we added a random mass scatter to each cluster, 
M' {a) = + 6\nM) by taking the fractional mass er¬ 

ror, cr(M^ a) ) / Mf ay for each cluster (see Table|lJ assuming the 
Gaussian distribution, and then computed the <^,_ sca]ing value 
for each realization. Here we used the mass errors for the gas 
mass proxy in Table |T] Note that the mean fractional error of 
50 clusters is about 0711, but here we included variations in 
the errors for different clusters. The figure shows that, if each 


cluster has a factor 2-3 larger mass error than the quoted error, 
the simulation results appear to reproduce the measurements. 
The mass errors of X-ray observables might underestimate the 
genuine mass uncertainty, perhaps due to the limitation of the 
X-ray based method or due to an unknown intrinsic scatter in 


the X-ray observable and halo mass relation (see also Okabe 
|et al.|2010b| for the similar discussion). 


3.4. The halo mass proxy relation of X-ray observables 


The method we have so far developed involves several as¬ 
sumptions. For instance, to implement the lensing stacking 
with NFW scaling, we need to assume several scaling rela¬ 
tions: the halo mass proxy relation of X-ray observables and 
the halo mass and concentration relation. In the following we 
address how possible variations in these scaling relations af¬ 
fect the NFW scaling results. 

First we study a possible bias in the X-ray inferred halo 
mass. Since the halo mass proxy relation of X-ray gas mass 
showed a better performance in the NFW scaling analysis (see 
Eq.|23]>, we here consider effects of possible variations in the 
gas mass relation on the results. To address this, we modify 
Eq. \22\ to parametrize the halo mass proxy relation as 


M, 


500 


10 14 Mc: 


= A x 11.6 x 


M. 


gas 


10 14 /r 3 / 2 M 0 


(24) 


where A is the normalization parameter and /3 is a slope pa¬ 
rameter of the halo mass dependence. The model with A — 1 
and /? = 1 corresponds to the self-similar model given by 
Eq. p2) . With varying the two parameters, A and /?, simul¬ 
taneously, we estimated the halo mass of each cluster from 
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Figure 16. Effects of variations in the halo mass proxy relation of X- 
ray gas mass on the NFW scaling results. We model the variations as 
M 5 ooc/[10 14 M q ] = A X 11.6 X (M gas /10 14 /!- 3/2 M Q y i > where A = 1 and 
P — 1 are our fiducial model corresponding to the self-similar scaling 
model. We estimated the best-fit parameters (the star symbol) by minimiz¬ 
ing the d 2 value with varying the normalization and mass slope parame¬ 
ters. The two contours correspond to the regions satisfying the conditions 
A d 1 = ^_ scaling (A >y 6) - ^*_ scaling W best ” flt ./8 best_flt ) = 2.3 or 6.17, respec- 
tively. The triangle symbol with errorbar denotes the result when varying the 
normalization parameter A alone, with fixing ft = 1. 


the X-ray gas mass based on the above proxy relation, and 
then redid the NFW scaling analysis. Figure fl6| shows the 
constraint regions in the two parameter space. TEe best-fit pa¬ 
rameters A = 0.86 ± 0.06 and /? = 0.66 ± 0.10, which has 
d 2 ,. = 493.8 compared to d 2 = 504.6 for the fidu- 

w-scaling r mm 


cial model as given by Eq. ( [23] ), corresponding to about 3cr 
improvement. Here we quoted the errorbars from the range 
Ad 2 < 1 with varying both A and / 3 , although the degeneracy 
between the two parameters is significant. Thus the scatters 
of 50 lensing distortion profiles prefer a weaker halo mass de¬ 
pendence than predicted by the self-similar scaling relation at 
a 3cr level. Note that, if the slope parameter is fixed to (3 — 1, 
i.e. the self-scaling relation, we obtained A = 1.05 ± 0.05 
(d 2 mn - 503.6, almost no change from the fiducical model). 
These results might be due to some residual uncertainty in 
our method, and would be worth further exploring by using a 
larger sample of clusters or an independent mass proxy rela¬ 
tion such as the Sunyaev-Zel’dovich effect. 


3.5. The halo mass and concentration relation 

Another important model ingredient in our analysis is the 
scaling relation of halo concentration with halo mass. We 
have so far employed the scaling relation in DK15 as for our 
default model. However, other works have proposed a dif¬ 
ferent scaling relation from DK15. For example, Duffy et al. 
(2008) proposed a different fitting formula of the c-M relation, 
and predicts a 20-30% lower concentration than in DK15 for 
cluster-scale halos: more exactly C 200 C ~ 3 or 4 for Duffy 
et al. ( j2008[ ) or DK15 for cluster-scale halos of several times 
10 14 /r'M G masses at z ~ 0 (see Figure 9 in DK15). How¬ 
ever, we found that, even if we use the scaling relation in 
Duffy et al. (2008) instead of DK15, it almost unchanges the 


d 2 ,■ value; more exactly, it enlarges the d 2 ,. value 

only by Ad 2 =* 1-2 for the HSE and gas mass proxy rela¬ 
tions. Hence the current data cannot discriminate these differ¬ 
ent models of c-M relation. 

Nevertheless one might ask whether or not our method al¬ 
lows us to constrain the underlying c-M relation. Assuming 
the parametrized form of c-M scaling relation given by 

J + P5) 

we minimized the d 2 ,. value with varying the normal- 
ization parameter and the mass slope parameter, f c and a. 
Here we took the halo mass inferred from the stacked lens¬ 
ing in Figure [6j M — 4 x 10 14 h 1 M G , for the pivot mass 
scale, and the redshift dependence is taken from |Duffy et al. 
(2008 Q Note that we fixed the mass normalization parameter 
to A = 1 for the halo mass proxy relation of X-ray observ¬ 
ables. Figure [17] shows that constraints on the two param¬ 
eters are significantly degenerate: the r^,_ scaling for the best- 
fit model is 526.3 or 503.5 for the HSE and gas mass, re¬ 
spectively, which is slightly smaller than our fiducial model, 
DK15, as found from Eq. (|23j ). The best-fit parameters are 
f c = 2.6^ and a = 0.08° ( ^7for HSE, while f c = 3.0 ± 0.3 
and a = -0.08 + 0.18 for the gas mass scaling relation. Thus 
the current data prefers the amplitude of concentration to be 
C 500 C — 2.6-3.0 for the 50 clusters of these mass scales, which 
is consistent with both the theory predictions in |Duffy et~ak 
( ]2008[ ) and DK15 within the errorbars, but cannot well con¬ 
strain the mass slope due to the limited statistics or a narrow 
range of halo masses. 


3.6. The halo mass profile 

The lensing analysis of NFW scaling rests on the assump¬ 
tion that the mass distribution in clusters follows the universal 
NFW profile. However, the NFW profile is the simplified pre¬ 
diction of A-body simulations, and a further improvement in 
our method might be available by employing a better model 
of the mass profile. 

Several works have pointed out variations in the inner re¬ 
gion of the mass profiles. For instance, there might be varia¬ 
tions in the inner slope of the mass profile (e.g., see Navarro 
et al. |2004[ and references therein). The baryonic processes 
would generally affect the inner structures, which tend to 
cause a greater mass concentration in the inner region and 
generally breaks the universality of the total mass profile. 
However, in this study, we looked into the cluster lensing sig¬ 
nals down to R ^ 0.14 h l Mpc and above, and these effects 
would be unlikely to largely change our results. 

Another interesting effect is a possible variation in the 
outer mass distribution at radii near to the virial radius or 
greater, as propos ed in |Diemer & Kravtsov| ( |2014| ) and |Ad-| 
hikari et ah (|2014J). These works claimed that the logarithmic 
slope of massive halos steepens more sharply than the NFW 
predicts, at the outer regions R > 0.5/Anoc, depending on the 
details of mass accretion and assembly history. This breaks 
to some extent the universality of NFW profile at these outer 
radii. Note that our analysis uses the lensing profile up to 
f?max = 2.8 Mpc/Vi, which corresponds to /? max ~ R 200 , for the 


7 Exactly speaking the fitting formula of Duffy et al. (2008) gives the c-M 
relation for M 2000 so we converted the scalingTelation to the relation between 
M 500 C and C 5000 and found that the redshift dependence is slightly modified 
from the original dependence (1 + z )~ 0 ' 47 by this conversion. 
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Figure 17. Effects of variations in the halo m ass a nd concentration relation on the NFW scaling results, for the HSE and gas mass cases, respectively. Here we 
parametrized the variations as c(M ) oc M a (Eq. |25| , and then minimized the d 2 -value with varying the normalization and mass slope parameters. The star symbol 
denotes the best-fit model. The two contours correspond to the regions satisfying the conditions Ad 2 = ^-scaling“ ^w-scaiing^ beSt_fit, ^ beSt_fit ) = ^-3 or 
6.17, respectively. The diamond symbol in each panel show s the parameters for our fiducial model DK15 at the mean redshift of clusters, z = 0.23, while the 
triangle symbol denotes the parameters of |Duffy et al.|(2008) . 


50 clusters. We tested the prediction of Diemer & Kravtsov 
( |2014| ) by using the fitting formula for a typical accretion his¬ 
tory that is kindly made available to us by Surhud More (jMore 
|et al.|2015) . However, we found that the current datasets can¬ 
not discriminate the steepened profile and the NFW profile at 
the outer radii. This would be interesting to further explore 
with an enlarged sample of clusters. 

3.7. A residual bias in source redshift 

As we discussed above, our results imply that the X-ray in¬ 
ferred mass may systematically underestimate the true mass: 
we found a possible bias of 5% level, although it is not sig¬ 
nificant (at a lcr level). Eq. ([7]) suggests that a 10% bias in 
halo mass corresponds to about 3% bias in the lensing ampli¬ 
tudes. This is a tiny amount, and may imply a residual error in 
the source redshift estimation. Due to the limited color infor¬ 
mation of the current data (mostly only 2 colors), we cannot 
resolve this, but a further study is definitely worth exploring. 
For the same reason, it is worth further looking into a possible 
remaining systematic error in the shape measurement. 


4. CONCLUSION AND DISCUSSION 

In this paper, we have developed a novel method of mea¬ 
suring the cluster lensing distortion profiles along the NFW 
prediction, one of the most important predictions of CDM 
structure formation model. The method measures the cluster 
lensing profiles by stacking the “scaled” amplitudes of back¬ 
ground galaxy ellipticities as a function of the “scaled” centric 
radius according to the NFW prediction of each cluster. To 
apply this method to real data, we combined the independent 
datasets for a nearly mass-selected sample of 50 massive clus¬ 
ters that are the Subaru weak lensing catalog in Okabe et al. 
(J2013) and their X-ray observables of XMM and/or Chandra 
satellites in Martino et al.](|2014|i. Here we used the X-ray ob¬ 
servables to infer the NFW parameters of each cluster; more 
precisely, we used the halo mass of each cluster based on the 


halo mass proxy relation of X-ray observables, either the hy¬ 
drostatic static equilibrium or the self-similar scaling relation 
of gas mass, and inferred the halo concentration from the c-M 
relation found in N-body simulations of DK15. We found a 4 
- 6<x level evidence of the existence of universal NFW profile 


23] and Figures [8] - [IT] ), 
ully studied a proper ra¬ 


in the 50 massive clusters (see Eq. 

To derive these results we have caref 
dial binning of the lensing distortion measurement and how to 
define the representative central value of each radial bin tak¬ 
ing into account the cluster-centric distances and the lensing 
weights of background galaxies in the annulus. Our results 
give a proof of concept of the method we developed in this 
paper. 

However, the improvement in the scatters of 50 cluster dis¬ 
tortion profiles due to the NFW scaling analysis is not as much 
as expected from theory using simulations of cluster based on 
high-resolution A-body simulations (Figures |~i~4|and [l5]l. We 
discussed that, in order to reconcile the difference between 
the measurements and the simulation expectation, we need to 
introduce additional halo mass scatters to each cluster, by an 
amount of cr(M)/M ~ 0.2-0.3 (see Figure [T5j). This implies 
intrinsi c scatt ers in the halo mass and X-ray observable rela¬ 
tion (jOkabe et al.[2010bjl. We also argued that the discrepancy 
might he due to an imper fect halo mass proxy relation of the 
X-ray observables (see § |3.4[ ). Hence it would be worth fur¬ 
ther exploring the method by combining different observables 
of clusters. A promising example is the Sunyaev-ZeTdovich 
(SZ) effect. By using or combining the X-ray, optical rich¬ 
ness and SZ effects to develop a well-calibrated relation be¬ 
tween halo mass and cluster observables for a suitable sample 
of massive clusters, we can explore a further improvement in 
constraining the universality of cluster mass distribution. In 
addition, we throughout used the model c-M relation to in¬ 
fer the halo concentration of each cluster. In other words, we 
ignored intrinsic scatters of halo concentration that is known 
to exist even for halos of a fixed mass scale from simulation 
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based studies. If we can use observables to estimate halo con¬ 
centration for each cluster, it might improve the NFW scaling 
results. For example, the concentration of member galaxies 
might be a good proxy of halo concentration on individual 
cluster basis. This would be worth exploring. 

Our method offers various applications. First, we inversely 
use the weak lensing analysis of NFW scaling to infer the un¬ 
derlying true relation between halo mass and cluster observ¬ 


ables. We made the initial attempt of this possibility in § 3.4 


Since the NFW scaling method up- or down-weights less or 
more massive clusters in order to make their profiles to be in 
the similar amplitudes, it can be applied to halos over a wider 
range of mass scales as long as the clusters in the sample 
follow the universal NFW profile. In this sense this method 
would be less sensitive to the selection effect of clusters in a 
sample. Secondly, we can use this method to explore the un¬ 
derlying true form of the halo mass profile or the ha lo m ass 
scaling relation with observables, as we attempted in §13.51 As 
claimed in|Diemer & Kravtsov|(]2014]), massive clusters might 
display a steeper profile at the outer radii around or beyond 
the virial radius than predicted by NFW model, depending on 
the mass accretion history. By subdividing clusters into sub¬ 
samples using a proxy to infer the mass accretion history, e.g. 
high or low concentration, we can use the NFW scaling analy¬ 
sis to explore the deviations from NFW prediction at the outer 
radii. This is a direct test of the hierarchical CDM structure 
formation model, and will be very interesting to explore. 

The weak lensing measurements of 50 massive clusters we 
used in this paper seem to be still limited by statistics, mainly 
due to a low number density of background galaxies, which 
we needed to take in order to define a secure sample of back¬ 
ground galaxies based on two passband data alone. Hence our 
method can be further improved by increasing background 
galaxies , based on photo-z inf ormation of more color infor¬ 
mation (Medezinski et al. 2013 ). We can also combine the 
lensing magnification bias measurement to improve the statis¬ 
tics. On-going wide-area optical surveys such as the HSC sur¬ 
vey and the DES survey promise to provide us with a much 
larger, well-calibrated sample of massive clusters, so it would 
be interesting to apply the method developed in this paper to 
those datasets in combination with other wavelength surveys 
such as X-ray or SZ effects. 
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